home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-19 / gpt32src.zip / STANDARD.C < prev    next >
C/C++ Source or Header  |  1992-03-25  |  17KB  |  807 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: standard.c,v 3.26 92/03/24 22:34:37 woo Exp Locker: woo $";
  3. #endif
  4.  
  5. /* GNUPLOT - standard.c */
  6. /*
  7.  * Copyright (C) 1986, 1987, 1990, 1991, 1992   Thomas Williams, Colin Kelley
  8.  *
  9.  * Permission to use, copy, and distribute this software and its
  10.  * documentation for any purpose with or without fee is hereby granted, 
  11.  * provided that the above copyright notice appear in all copies and 
  12.  * that both that copyright notice and this permission notice appear 
  13.  * in supporting documentation.
  14.  *
  15.  * Permission to modify the software is granted, but not the right to
  16.  * distribute the modified code.  Modifications are to be distributed 
  17.  * as patches to released version.
  18.  *  
  19.  * This software is provided "as is" without express or implied warranty.
  20.  * 
  21.  *
  22.  * AUTHORS
  23.  * 
  24.  *   Original Software:
  25.  *     Thomas Williams,  Colin Kelley.
  26.  * 
  27.  *   Gnuplot 2.0 additions:
  28.  *       Russell Lang, Dave Kotz, John Campbell.
  29.  *
  30.  *   Gnuplot 3.0 additions:
  31.  *       Gershon Elber and many others.
  32.  * 
  33.  * Send your comments or suggestions to 
  34.  *  info-gnuplot@ames.arc.nasa.gov.
  35.  * This is a mailing list; to join it send a note to 
  36.  *  info-gnuplot-request@ames.arc.nasa.gov.  
  37.  * Send bug reports to
  38.  *  bug-gnuplot@ames.arc.nasa.gov.
  39.  */
  40.  
  41. #include <math.h>
  42. #include <stdio.h>
  43. #include "plot.h"
  44.  
  45. #ifdef vms
  46. #include <errno.h>
  47. #else
  48. extern int errno;
  49. #endif /* vms */
  50.  
  51.  
  52. extern struct value stack[STACK_DEPTH];
  53. extern int s_p;
  54. extern double zero;
  55.  
  56. struct value *pop(), *complex(), *integer();
  57.  
  58. double magnitude(), angle(), real(), imag();
  59.  
  60. /* The bessel function approximations here are from
  61.  * "Computer Approximations"
  62.  * by Hart, Cheney et al.
  63.  * John Wiley & Sons, 1968
  64.  */
  65.  
  66. /* There appears to be a mistake in Hart, Cheney et al. on page 149.
  67.  * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
  68.  *               Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
  69.  * In the functions below, Qn(x) is implementated using the later
  70.  * equation.
  71.  * These bessel functions are accurate to about 1e-13
  72.  */
  73.  
  74. #define PI_ON_FOUR       0.78539816339744830961566084581987572
  75. #define PI_ON_TWO        1.57079632679489661923131269163975144
  76. #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
  77. #define TWO_ON_PI        0.63661977236758134307553505349005744
  78.  
  79. static double dzero = 0.0;
  80.  
  81. /* jzero for x in [0,8]
  82.  * Index 5849, 19.22 digits precision
  83.  */
  84. static double pjzero[] = {
  85.      0.4933787251794133561816813446e+21,
  86.     -0.11791576291076105360384408e+21,
  87.      0.6382059341072356562289432465e+19,
  88.     -0.1367620353088171386865416609e+18,
  89.      0.1434354939140346111664316553e+16,
  90.     -0.8085222034853793871199468171e+13,
  91.      0.2507158285536881945555156435e+11,
  92.     -0.4050412371833132706360663322e+8,
  93.      0.2685786856980014981415848441e+5
  94. };
  95.  
  96. static double qjzero[] = {
  97.     0.4933787251794133562113278438e+21,
  98.     0.5428918384092285160200195092e+19,
  99.     0.3024635616709462698627330784e+17,
  100.     0.1127756739679798507056031594e+15,
  101.     0.3123043114941213172572469442e+12,
  102.     0.669998767298223967181402866e+9,
  103.     0.1114636098462985378182402543e+7,
  104.     0.1363063652328970604442810507e+4,
  105.     0.1e+1
  106. };
  107.  
  108. /* pzero for x in [8,inf]
  109.  * Index 6548, 18.16 digits precision
  110.  */
  111. static double ppzero[] = {
  112.     0.2277909019730468430227002627e+5,
  113.     0.4134538663958076579678016384e+5,
  114.     0.2117052338086494432193395727e+5,
  115.     0.348064864432492703474453111e+4,
  116.     0.15376201909008354295771715e+3,
  117.     0.889615484242104552360748e+0
  118. };
  119.  
  120. static double qpzero[] = {
  121.     0.2277909019730468431768423768e+5,
  122.     0.4137041249551041663989198384e+5,
  123.     0.2121535056188011573042256764e+5,
  124.     0.350287351382356082073561423e+4,
  125.     0.15711159858080893649068482e+3,
  126.     0.1e+1
  127. };
  128.  
  129. /* qzero for x in [8,inf]
  130.  * Index 6948, 18.33 digits precision
  131.  */
  132. static double pqzero[] = {
  133.     -0.8922660020080009409846916e+2,
  134.     -0.18591953644342993800252169e+3,
  135.     -0.11183429920482737611262123e+3,
  136.     -0.2230026166621419847169915e+2,
  137.     -0.124410267458356384591379e+1,
  138.     -0.8803330304868075181663e-2,
  139. };
  140.  
  141. static double qqzero[] = {
  142.     0.571050241285120619052476459e+4,
  143.     0.1195113154343461364695265329e+5,
  144.     0.726427801692110188369134506e+4,
  145.     0.148872312322837565816134698e+4,
  146.     0.9059376959499312585881878e+2,
  147.     0.1e+1
  148. };
  149.  
  150.  
  151. /* yzero for x in [0,8]
  152.  * Index 6245, 18.78 digits precision
  153.  */
  154. static double pyzero[] = {
  155.     -0.2750286678629109583701933175e+20,
  156.      0.6587473275719554925999402049e+20,
  157.     -0.5247065581112764941297350814e+19,
  158.      0.1375624316399344078571335453e+18,
  159.     -0.1648605817185729473122082537e+16,
  160.      0.1025520859686394284509167421e+14,
  161.     -0.3436371222979040378171030138e+11,
  162.      0.5915213465686889654273830069e+8,
  163.     -0.4137035497933148554125235152e+5
  164. };
  165.  
  166. static double qyzero[] = {
  167.     0.3726458838986165881989980739e+21,
  168.     0.4192417043410839973904769661e+19,
  169.     0.2392883043499781857439356652e+17,
  170.     0.9162038034075185262489147968e+14,
  171.     0.2613065755041081249568482092e+12,
  172.     0.5795122640700729537380087915e+9,
  173.     0.1001702641288906265666651753e+7,
  174.     0.1282452772478993804176329391e+4,
  175.     0.1e+1
  176. };
  177.  
  178.  
  179. /* jone for x in [0,8]
  180.  * Index 6050, 20.98 digits precision
  181.  */
  182. static double pjone[] = {
  183.      0.581199354001606143928050809e+21,
  184.     -0.6672106568924916298020941484e+20,
  185.      0.2316433580634002297931815435e+19,
  186.     -0.3588817569910106050743641413e+17,
  187.      0.2908795263834775409737601689e+15,
  188.     -0.1322983480332126453125473247e+13,
  189.      0.3413234182301700539091292655e+10,
  190.     -0.4695753530642995859767162166e+7,
  191.      0.270112271089232341485679099e+4
  192. };
  193.  
  194. static double qjone[] = {
  195.     0.11623987080032122878585294e+22,
  196.     0.1185770712190320999837113348e+20,
  197.     0.6092061398917521746105196863e+17,
  198.     0.2081661221307607351240184229e+15,
  199.     0.5243710262167649715406728642e+12,
  200.     0.1013863514358673989967045588e+10,
  201.     0.1501793594998585505921097578e+7,
  202.     0.1606931573481487801970916749e+4,
  203.     0.1e+1
  204. };
  205.  
  206.  
  207. /* pone for x in [8,inf]
  208.  * Index 6749, 18.11 digits precision
  209.  */
  210. static double ppone[] = {
  211.     0.352246649133679798341724373e+5,
  212.     0.62758845247161281269005675e+5,
  213.     0.313539631109159574238669888e+5,
  214.     0.49854832060594338434500455e+4,
  215.     0.2111529182853962382105718e+3,
  216.     0.12571716929145341558495e+1
  217. };
  218.  
  219. static double qpone[] = {
  220.     0.352246649133679798068390431e+5,
  221.     0.626943469593560511888833731e+5,
  222.     0.312404063819041039923015703e+5,
  223.     0.4930396490181088979386097e+4,
  224.     0.2030775189134759322293574e+3,
  225.     0.1e+1
  226. };
  227.  
  228. /* qone for x in [8,inf]
  229.  * Index 7149, 18.28 digits precision
  230.  */
  231. static double pqone[] = {
  232.     0.3511751914303552822533318e+3,
  233.     0.7210391804904475039280863e+3,
  234.     0.4259873011654442389886993e+3,
  235.     0.831898957673850827325226e+2,
  236.     0.45681716295512267064405e+1,
  237.     0.3532840052740123642735e-1
  238. };
  239.  
  240. static double qqone[] = {
  241.     0.74917374171809127714519505e+4,
  242.     0.154141773392650970499848051e+5,
  243.     0.91522317015169922705904727e+4,
  244.     0.18111867005523513506724158e+4,
  245.     0.1038187585462133728776636e+3,
  246.     0.1e+1
  247. };
  248.  
  249.  
  250. /* yone for x in [0,8]
  251.  * Index 6444, 18.24 digits precision
  252.  */
  253. static double pyone[] = {
  254.     -0.2923821961532962543101048748e+20,
  255.      0.7748520682186839645088094202e+19,
  256.     -0.3441048063084114446185461344e+18,
  257.      0.5915160760490070618496315281e+16,
  258.     -0.4863316942567175074828129117e+14,
  259.      0.2049696673745662182619800495e+12,
  260.     -0.4289471968855248801821819588e+9,
  261.      0.3556924009830526056691325215e+6
  262. };
  263.  
  264. static double qyone[] = {
  265.     0.1491311511302920350174081355e+21,
  266.     0.1818662841706134986885065935e+19,
  267.     0.113163938269888452690508283e+17,
  268.     0.4755173588888137713092774006e+14,
  269.     0.1500221699156708987166369115e+12,
  270.     0.3716660798621930285596927703e+9,
  271.     0.726914730719888456980191315e+6,
  272.     0.10726961437789255233221267e+4,
  273.     0.1e+1
  274. };
  275.  
  276.  
  277. f_real()
  278. {
  279. struct value a;
  280.     push( complex(&a,real(pop(&a)), 0.0) );
  281. }
  282.  
  283. f_imag()
  284. {
  285. struct value a;
  286.     push( complex(&a,imag(pop(&a)), 0.0) );
  287. }
  288.  
  289. f_arg()
  290. {
  291. struct value a;
  292.     push( complex(&a,angle(pop(&a)), 0.0) );
  293. }
  294.  
  295. f_conjg()
  296. {
  297. struct value a;
  298.     (void) pop(&a);
  299.     push( complex(&a,real(&a),-imag(&a) ));
  300. }
  301.  
  302. f_sin()
  303. {
  304. struct value a;
  305.     (void) pop(&a);
  306.     push( complex(&a,sin(real(&a))*cosh(imag(&a)), cos(real(&a))*sinh(imag(&a))) );
  307. }
  308.  
  309. f_cos()
  310. {
  311. struct value a;
  312.     (void) pop(&a);
  313.     push( complex(&a,cos(real(&a))*cosh(imag(&a)), -sin(real(&a))*sinh(imag(&a))));
  314. }
  315.  
  316. f_tan()
  317. {
  318. struct value a;
  319. register double den;
  320.     (void) pop(&a);
  321.     if (imag(&a) == 0.0)
  322.         push( complex(&a,tan(real(&a)),0.0) );
  323.     else {
  324.         den = cos(2*real(&a))+cosh(2*imag(&a));
  325.         if (den == 0.0) {
  326.             undefined = TRUE;
  327.             push( &a );
  328.         }
  329.         else
  330.             push( complex(&a,sin(2*real(&a))/den, sinh(2*imag(&a))/den) );
  331.     }
  332. }
  333.  
  334. f_asin()
  335. {
  336. struct value a;
  337. register double alpha, beta, x, y;
  338.     (void) pop(&a);
  339.     x = real(&a); y = imag(&a);
  340.     if (y == 0.0) {
  341.         if (fabs(x) > 1.0) {
  342.             undefined = TRUE;
  343.             push(complex(&a,0.0, 0.0));
  344.         } else
  345.             push( complex(&a,asin(x),0.0) );
  346.     } else {
  347.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  348.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  349.         push( complex(&a,asin(beta), log(alpha + sqrt(alpha*alpha-1))) );
  350.     }
  351. }
  352.  
  353. f_acos()
  354. {
  355. struct value a;
  356. register double alpha, beta, x, y;
  357.     (void) pop(&a);
  358.     x = real(&a); y = imag(&a);
  359.     if (y == 0.0) {
  360.         if (fabs(x) > 1.0) {
  361.             undefined = TRUE;
  362.             push(complex(&a,0.0, 0.0));
  363.         } else
  364.             push( complex(&a,acos(x),0.0) );
  365.     } else {
  366.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  367.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  368.         push( complex(&a,acos(beta), log(alpha + sqrt(alpha*alpha-1))) );
  369.     }
  370. }
  371.  
  372. f_atan()
  373. {
  374. struct value a;
  375. register double x, y, u, v, w, z;
  376.     (void) pop(&a);
  377.     x = real(&a); y = imag(&a);
  378.     if (y == 0.0)
  379.         push( complex(&a,atan(x), 0.0) );
  380.     else if (x == 0.0 && fabs(y) == 1.0) {
  381.         undefined = TRUE;
  382.         push(complex(&a,0.0, 0.0));
  383.     } else {
  384.             if (x >= 0) {
  385.                 u = x;
  386.             v = y;
  387.         } else {
  388.                 u = -x;
  389.             v = -y;
  390.         }
  391.         
  392.             z = atan(2*u/(1-u*u-v*v));
  393.         w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
  394.         if (z < 0)
  395.                 z = z + 2*PI_ON_TWO;
  396.         if (x < 0) {
  397.                 z = -z;
  398.             w = -w;
  399.         }
  400.         push( complex(&a,0.5*z, w) );
  401.     }
  402. }
  403.  
  404. f_sinh()
  405. {
  406. struct value a;
  407.     (void) pop(&a);
  408.     push( complex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
  409. }
  410.  
  411. f_cosh()
  412. {
  413. struct value a;
  414.     (void) pop(&a);
  415.     push( complex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
  416. }
  417.  
  418. f_tanh()
  419. {
  420. struct value a;
  421. register double den;
  422.     (void) pop(&a);
  423.     den = cosh(2*real(&a)) + cos(2*imag(&a));
  424.     push( complex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
  425. }
  426.  
  427. f_int()
  428. {
  429. struct value a;
  430.     push( integer(&a,(int)real(pop(&a))) );
  431. }
  432.  
  433.  
  434. f_abs()
  435. {
  436. struct value a;
  437.     (void) pop(&a);
  438.     switch (a.type) {
  439.         case INT:
  440.             push( integer(&a,abs(a.v.int_val)) );            
  441.             break;
  442.         case CMPLX:
  443.             push( complex(&a,magnitude(&a), 0.0) );
  444.     }
  445. }
  446.  
  447. f_sgn()
  448. {
  449. struct value a;
  450.     (void) pop(&a);
  451.     switch(a.type) {
  452.         case INT:
  453.             push( integer(&a,(a.v.int_val > 0) ? 1 : 
  454.                     (a.v.int_val < 0) ? -1 : 0) );
  455.             break;
  456.         case CMPLX:
  457.             push( integer(&a,(a.v.cmplx_val.real > 0.0) ? 1 : 
  458.                     (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
  459.             break;
  460.     }
  461. }
  462.  
  463.  
  464. f_sqrt()
  465. {
  466. struct value a;
  467. register double mag, ang;
  468.     (void) pop(&a);
  469.     mag = sqrt(magnitude(&a));
  470.     if (imag(&a) == 0.0 && real(&a) < 0.0)
  471.         push( complex(&a,0.0,mag) );
  472.     else
  473.     {
  474.         if ( (ang = angle(&a)) < 0.0)
  475.             ang += 2*Pi;
  476.         ang /= 2;
  477.         push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  478.     }
  479. }
  480.  
  481.  
  482. f_exp()
  483. {
  484. struct value a;
  485. register double mag, ang;
  486.     (void) pop(&a);
  487.     mag = exp(real(&a));
  488.     ang = imag(&a);
  489.     push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  490. }
  491.  
  492.  
  493. f_log10()
  494. {
  495. struct value a;
  496. register double l10;;
  497.     (void) pop(&a);
  498.     l10 = log(10.0);    /***** replace with a constant! ******/
  499.     push( complex(&a,log(magnitude(&a))/l10, angle(&a)/l10) );
  500. }
  501.  
  502.  
  503. f_log()
  504. {
  505. struct value a;
  506.     (void) pop(&a);
  507.     push( complex(&a,log(magnitude(&a)), angle(&a)) );
  508. }
  509.  
  510.  
  511. f_floor()
  512. {
  513. struct value a;
  514.  
  515.     (void) pop(&a);
  516.     switch (a.type) {
  517.         case INT:
  518.             push( integer(&a,(int)floor((double)a.v.int_val)));            
  519.             break;
  520.         case CMPLX:
  521.             push( integer(&a,(int)floor(a.v.cmplx_val.real)));
  522.     }
  523. }
  524.  
  525.  
  526. f_ceil()
  527. {
  528. struct value a;
  529.  
  530.     (void) pop(&a);
  531.     switch (a.type) {
  532.         case INT:
  533.             push( integer(&a,(int)ceil((double)a.v.int_val)));            
  534.             break;
  535.         case CMPLX:
  536.             push( integer(&a,(int)ceil(a.v.cmplx_val.real)));
  537.     }
  538. }
  539.  
  540. #ifdef GAMMA
  541.  
  542. f_gamma()
  543. {
  544. extern int signgam;
  545. register double y;
  546. struct value a;
  547.  
  548.     y = GAMMA(real(pop(&a)));
  549.     if (y > 88.0) {
  550.         undefined = TRUE;
  551.         push( integer(&a,0) );
  552.     }
  553.     else
  554.         push( complex(&a,signgam * exp(y),0.0) );
  555. }
  556.  
  557. #endif /* GAMMA */
  558.  
  559.  
  560. /* bessel function approximations */
  561. double jzero(x)
  562. double x;
  563. {
  564. double p, q, x2;
  565. int n;
  566.  
  567.     x2 = x * x;
  568.     p = pjzero[8];
  569.     q = qjzero[8];
  570.     for (n=7; n>=0; n--) {
  571.         p = p*x2 + pjzero[n];
  572.         q = q*x2 + qjzero[n];
  573.     }
  574.     return(p/q);
  575. }
  576.  
  577. double pzero(x)
  578. double x;
  579. {
  580. double p, q, z, z2;
  581. int n;
  582.  
  583.     z = 8.0 / x;
  584.     z2 = z * z;
  585.     p = ppzero[5];
  586.     q = qpzero[5];
  587.     for (n=4; n>=0; n--) {
  588.         p = p*z2 + ppzero[n];
  589.         q = q*z2 + qpzero[n];
  590.     }
  591.     return(p/q);
  592. }
  593.  
  594. double qzero(x)
  595. double x;
  596. {
  597. double p, q, z, z2;
  598. int n;
  599.  
  600.     z = 8.0 / x;
  601.     z2 = z * z;
  602.     p = pqzero[5];
  603.     q = qqzero[5];
  604.     for (n=4; n>=0; n--) {
  605.         p = p*z2 + pqzero[n];
  606.         q = q*z2 + qqzero[n];
  607.     }
  608.     return(p/q);
  609. }
  610.  
  611. double yzero(x)
  612. double x;
  613. {
  614. double p, q, x2;
  615. int n;
  616.  
  617.     x2 = x * x;
  618.     p = pyzero[8];
  619.     q = qyzero[8];
  620.     for (n=7; n>=0; n--) {
  621.         p = p*x2 + pyzero[n];
  622.         q = q*x2 + qyzero[n];
  623.     }
  624.     return(p/q);
  625. }
  626.  
  627. double rj0(x)
  628. double x;
  629. {
  630.     if ( x <= 0.0 )
  631.         x = -x;
  632.     if ( x < 8.0 )
  633.         return(jzero(x));
  634.     else
  635.         return( sqrt(TWO_ON_PI/x) *
  636.             (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
  637.  
  638. }
  639.  
  640. double ry0(x)
  641. double x;
  642. {
  643.     if ( x < 0.0 )
  644.         return(dzero/dzero); /* error */
  645.     if ( x < 8.0 )
  646.         return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
  647.     else
  648.         return( sqrt(TWO_ON_PI/x) *
  649.             (pzero(x)*sin(x-PI_ON_FOUR) + 
  650.             (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
  651.  
  652. }
  653.  
  654.  
  655. double jone(x)
  656. double x;
  657. {
  658. double p, q, x2;
  659. int n;
  660.  
  661.     x2 = x * x;
  662.     p = pjone[8];
  663.     q = qjone[8];
  664.     for (n=7; n>=0; n--) {
  665.         p = p*x2 + pjone[n];
  666.         q = q*x2 + qjone[n];
  667.     }
  668.     return(p/q);
  669. }
  670.  
  671. double pone(x)
  672. double x;
  673. {
  674. double p, q, z, z2;
  675. int n;
  676.  
  677.     z = 8.0 / x;
  678.     z2 = z * z;
  679.     p = ppone[5];
  680.     q = qpone[5];
  681.     for (n=4; n>=0; n--) {
  682.         p = p*z2 + ppone[n];
  683.         q = q*z2 + qpone[n];
  684.     }
  685.     return(p/q);
  686. }
  687.  
  688. double qone(x)
  689. double x;
  690. {
  691. double p, q, z, z2;
  692. int n;
  693.  
  694.     z = 8.0 / x;
  695.     z2 = z * z;
  696.     p = pqone[5];
  697.     q = qqone[5];
  698.     for (n=4; n>=0; n--) {
  699.         p = p*z2 + pqone[n];
  700.         q = q*z2 + qqone[n];
  701.     }
  702.     return(p/q);
  703. }
  704.  
  705. double yone(x)
  706. double x;
  707. {
  708. double p, q, x2;
  709. int n;
  710.  
  711.     x2 = x * x;
  712.     p = 0.0;
  713.     q = qyone[8];
  714.     for (n=7; n>=0; n--) {
  715.         p = p*x2 + pyone[n];
  716.         q = q*x2 + qyone[n];
  717.     }
  718.     return(p/q);
  719. }
  720.  
  721. double rj1(x)
  722. double x;
  723. {
  724. double v,w;
  725.     v = x;
  726.     if ( x < 0.0 )
  727.         x = -x;
  728.     if ( x < 8.0 )
  729.         return(v*jone(x));
  730.     else {
  731.         w = sqrt(TWO_ON_PI/x) *
  732.             (pone(x)*cos(x-THREE_PI_ON_FOUR) - 
  733.                8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
  734.         if (v < 0.0)
  735.             w = -w;
  736.         return( w );
  737.     }
  738. }
  739.  
  740. double ry1(x)
  741. double x;
  742. {
  743.     if ( x <= 0.0 )
  744.         return(dzero/dzero); /* error */
  745.     if ( x < 8.0 )
  746.         return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
  747.     else
  748.         return( sqrt(TWO_ON_PI/x) *
  749.             (pone(x)*sin(x-THREE_PI_ON_FOUR) + 
  750.             (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
  751. }
  752.  
  753.  
  754. f_besj0()    
  755. {
  756. struct value a;
  757. double x;
  758.     (void) pop(&a);
  759.     if (imag(&a) > zero)
  760.         int_error("can only do bessel functions of reals",NO_CARET);
  761.     push( complex(&a,rj0(real(&a)),0.0) );
  762. }
  763.  
  764.  
  765. f_besj1()    
  766. {
  767. struct value a;
  768. double x;
  769.     (void) pop(&a);
  770.     if (imag(&a) > zero)
  771.         int_error("can only do bessel functions of reals",NO_CARET);
  772.     push( complex(&a,rj1(real(&a)),0.0) );
  773. }
  774.  
  775.  
  776. f_besy0()    
  777. {
  778. struct value a;
  779. double x;
  780.     (void) pop(&a);
  781.     if (imag(&a) > zero)
  782.         int_error("can only do bessel functions of reals",NO_CARET);
  783.     if (real(&a) > 0.0)
  784.         push( complex(&a,ry0(real(&a)),0.0) );
  785.     else {
  786.         push( complex(&a,0.0,0.0) );
  787.         undefined = TRUE ;
  788.     }
  789. }
  790.  
  791.  
  792. f_besy1()    
  793. {
  794. struct value a;
  795. double x;
  796.     (void) pop(&a);
  797.     if (imag(&a) > zero)
  798.         int_error("can only do bessel functions of reals",NO_CARET);
  799.     if (real(&a) > 0.0)
  800.         push( complex(&a,ry1(real(&a)),0.0) );
  801.     else {
  802.         push( complex(&a,0.0,0.0) );
  803.         undefined = TRUE ;
  804.     }
  805. }
  806.  
  807.